A multimode model for projective photon-counting measurements 



Rosa Tualle-Brouri ,* Alexei Ourjoumtsev, Aurelien Dantan, and Philippe Grangier 
Laboratoire Charles Fabry de I'Institut d'Optique, C'NRS UMR 8501, 
Universite Paris Sud XI, 91127 Palaiseau, France 



0^ 
O 
O 
(N 



> 

00 
(N 

(N 
O 

o 



X 



Martijn Wubs^'^ and Anders S. S0rensen§ 
t Niels Bohr Inter-national Academy & § QUANTOP 
The Niels Bohr Institute, Blegdamsvej 17, DK-2100 Copenhagen, Denmark 
(Dated: Version 17, date: February 17, 2009) 

We present a general model to account for the multimode nature of the quantum electromagnetic 
field in projective photon-counting measurements. We focus on photon-subtraction experiments, 
where non-gaussian states are produced conditionally. These are useful states for continuous- variable 
quantum information processing. We present a general method called mode reduction that reduces 
the multimode model to an effective two-mode problem. We apply this method to a multimode 
model describing broadband parametric downconversion, thereby improving the analysis of existing 
experimental results. The main improvement is that spatial and frequency filters before the photon 
detector are taken into account explicitly. We find excellent agreement with previously published 
experimental results, using fewer free parameters than before, and discuss the implications of our 
analysis for the optimized production of states with negative Wigner functions. 

PACS numbers: 03.67.-a, 42.50.Dv, 03.65.Wj 



I. INTRODUCTION 

The ability to prepare and measure specific quantum 
states of the light is the keystone of many quantum infor- 
mation processing (QIP) protocols. These states can be 
described either with discrete variables in terms of pho- 
tons, or with continuous variables in terms of waves. In 
the latter case, the physical quantities of interest are the 
amplitude and the phase of the light wave, or their Carte- 
sian counterparts called quadratures x and p. A very 
convenient representation of the quantum state is then 
provided by the Wigner function W{x,p)^ which corre- 
sponds to a quasi-probability distribution of the quadra- 
tures, 'quasi-' because W may assume negative values. 

An important task for QIP is the ability to undo effects 
of decoherence by 'distillation': to obtain a single quan- 
tum state that is more pure from two or more copies 
that have undergone decoherence. Since states of light 
with gaussian Wigner functions cannot be distilled with 
gaussian operations [1, 2], one is left with two strate- 
gies: either to distill gaussian states with non-gaussian 
operations, or to distill non-gaussian states with gaus- 
sian operations [3]. This paper is a contribution to the 
latter strategy, and focuses on the preparation of the non- 
gaussian states rather than on their distillation. 

The negativity of the Wigner function is a standard 
figure of merit, quantifying at the same time how non- 
gaussian and how non-classical a quantum state is [4, 5] . 
One way of obtaining non-gausian states is by condi- 
tional photon-counting measurements, as first proposed 
by Dakna et al. [6] . It was soon realized that such condi- 
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tional measurements can improve quantum teleportation 
of continuous variables [7]. In recent years, several exper- 
iments [8, 9, 10, 11, 12, 13, 14, 15] combining continuous- 
and discrete-variable tools allowed for preparing and ob- 
serving quantum states of free-propagating light with 
negative Wigner functions [16]. 

Many of these experiments are based on the use of a 
squeezed vacuum produced by parametric fluorescence, 
which involves many optical modes [17]. This multimode 
nature is exhibited in both continuous-wave (CW) op- 
eration, using optical parametric oscillators (OPO) be- 
low threshold [18, 19], and in pulsed experiments with a 
single-pass high amplification. In order to make accurate 
predictions, it is therefore crucial to develop multimode 
theoretical models. This was done in [18, 19, 20, 21] for 
setups using an OPO, and in [22, 23] for pulsed 1-photon 
Fock state tomography in a case of very low squeezing, 
when Fock states expansion are limited to 1 photon only. 
However, these models do not fully account for all phe- 
nomena linked to the non-constant space and time pro- 
files of the modes under study, the spatial pulse profile in 
the transverse direction for example, and they especially 
do not account for gain-induced distortions in the para- 
metric amplification process [24] . As we will see later on, 
these phenomena are one main signature of this multi- 
mode nature, and are critical in the case of single-pass 
pulsed experiments. 

In this paper we propose an alternative general frame- 
work to describe the generation of squeezed light, with a 
twofold goal: first, to introduce a method that reduces a 
multimode model to an effective two-mode description. 
Second, to show that a specific spatio-temporal multi- 
mode model for photon-subtraction experiments and our 
mode-reduction analysis thereof give an improved un- 
derstanding of state-of-the-art photon-subtraction exper- 
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FIG. 1: Projective photon-counting measurement; in our 
model we assume that we start with a multimode vacuum 
input state which is subject to a general multimode Bogoli- 
ubov transform U. After the transformation a few modes are 
filtered out and result in a detection event in the avalanche 
photodiode (APD). Such a detection event in the APD pre- 
pares a state in the single mode which is mode analyzed by 
homodyne detection. 



iments. 

In Sec. II we show how to reduce a multi-mode model 

to an cfFcctivc two-mode model. This mode-reduction 
procedure is then applied in Sec. Ill to give an im- 
proved analysis of the photon-subtraction experiments 
of Ref. [10]. We discuss the method and its application 
and conclude in Sec. IV. Some technicalities are deferred 
to two Appendices. 



II. REDUCTION OF MULTIMODE MODEL 

A. General multimode model 

As a starting point we have a complete set of spatial, 
temporal, or spatio-temporal optical modes, in terms of 
which the light propagation can be described. The modes 
have field operators aj and aj with bosonic commutation 
relations f 

fe] ~ ^jk- The a,j may also stand for contin- 
uous operators like a(t), where t is time, with commuta- 
tion relation [a(i), a^(i')] = 6{t — t'), in which case sums 
over modes are to be replaced with integrals. To simplify 
the notation even further, we introduce a, the column 
vector of all the aj . Similarly a^ is the row vector which 
has operators aj as components. We shall also need the 

row vector a^ and the column vector a* = a'^'^. Similar 
notation will be used for other vectors. 

In this section we consider a general unitary transfor- 
mation U in which output quadratures linearly depend on 
input quadratures, thus preserving the gaussian nature of 
the quantum fields. In fact, the only non-gaussian oper- 
ation will be the projective measurement, corresponding 
to a detection event in a subset of the output modes by 
an avalanche photodiode (APD), see Figure 1. We will 
assume the evolution operator U in Fig. 1 to be a Bogoli- 
ubov transformation, where the output field operators 



depend linearly on the input ones: 

aout = all = u a + V a* , 
where u and v are two matrices which satisfy 



UW — w 
uv^ — vu^ = 



u'u — v'^v = 1, 



(1) 



(2a) 
(2b) 



whereby the commutation relations of the field opera- 
tors are preserved. The output state of the light can 
be characterized by doing homodyne measurements on a 
normalized mode described by tphj which has the mode 
operator 



ah = tfi a. 



(3) 



This mode can be defined as the mode that perfectly 
matches the local oscillator of the homodyne detector. 
After the Bogoliubov transform (1) it is described by 



ah,out = U'^ahU = iljl{ua + va*). 



(4) 



Without conditioning upon photon detection, and since 
the gaussian nature of the initial vacuum state is pre- 
served by the Bogoliubov transform, the homodyne mea- 
surements will show gaussian Wigner functions corre- 
sponding to squeezed vacimm states of light. 

The point is now how to describe the output state con- 
ditional upon detection of a photon by the APD, which is 
a projective measurement. Analogous to the homodyne 
detection, we assume the photon detection mode j to be 
described by a normalized state (pd{j), with correspond- 
ing field operator 



(5) 



After the time evolution described by Eq. (1), the output 
at the photon detector becomes: 



ad,out(i) = U^ad{j)U = ^^{j){ua + va*) 



(6) 



If we only know that a photon has been detected but 

not in which detection mode, then we should average the 
conditional output state over all these detection modes, 
as is shown in more detail below. 

At this point we have a multimode output state that in 
principle we can update given the detection of a photon 
in the detector j. In practice it is convenient to first 
simplify the multimode expressions (4) and (6). 



B. Mode reduction 

The first and central step in the mode-reduction pro- 
cedure is to rewrite the homodyne mode Eq. (4) in the 
form 



a;i,out = r? ao + a 



■pal 



(7) 
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FIG. 2: Equivalent model for the Bogoliubov transform U: 
a perfect single-mode degenerate optical parametric amplifier 

(DOPA) with squeezing parameter r is followed by a perfect 
non-degenerate optical parametric amplifier (NDOPA) with 
squeezing parameter gr. 



in terms of new mode operators Cq^J with standard com- 
mutation relations. The coefficients r], a, and /3 in Eq. (7) 
are found as follows. Besides having standard commuta- 
tions, So and ai in Eq. (7) should annihilate the vacmim 
state, so that ao must contain all annihilation operators 
in Eq. (4): 



T]ao = tpj^u a. 



(8) 



with 77 fixed up to a phase factor by [ao,aJ] = 1. We 
choose T] to be real- valued and positive so that 



Furthermore, from [ao,a|] =0 it follows that 



(9) 



(10) 



A complex value for a can be removed with a redefinition 
of the phase of the homodyne mode, which just means 

that we can assume that the state is squeezed in the x or 
p direction. Finally, /3 can be found from the fact that 
ah,out and dj^ have commutator one, as have ai and 



(11) 



Here again we used the freedom to choose (3 real- valued 
and nonnegative. This completes the mathematics of the 
mode reduction of the multimode homodyning signal. 

The physical argument that two and only two modes 
should remain goes as follows. The squeezed vacuum af- 
ter the Bogoliubov transform can only be a centered gaus- 
sian state, hence it is fully described by only the variances 
Vx and Vp. The squeezed vacuum output can therefore 
be modeled [25, 26] by a perfect single-mode degenerate 
optical parametric amplifier (DOPA) with squeezing pa- 
rameter r, followed by a perfect non-degenerate optical 
parametric amplifier (NDOPA) with squeezing parame- 
ter gr, as presented in Figure 2. This is a two-mode 



model, with a Hilbert space 7^2- The values of the param- 
eters r and g can be deduced from the two independent 
coefficients in Eqs. (9-11), using: 



T] = cosh(r) cosh(£fr) 
a = sinh(r) cosh(£fr) 
(3 = smh{gr). 



(12a) 
(12b) 
(12c) 



A related squeezing parameter that we will also use in 
the following is s = exp(— 2r). 



C. Conditioning upon photon detection 

Wc now condition upon the measurement of a click in 
the photon detector (APD). We assume to be in the limit 
that the average number of photons per pulse entering 
the photon detector is much less than one. Then a single 
click in the detector corresponds to the detection of a 
single photon. 

One can make a reduced-mode description of the pho- 
ton detection operator ad.out{j) of Eq. (5), analogous to 
Eq. (7). The operator (id.outij) can be expanded into a 
part acting on H2, plus a component ad±{j) acting on 
the complementary space 'H± orthogonal to W2: 

ad,out{j) = Ijal + Sja\ + ejao + KjCLi + ad±{j). (13) 

Note that a4±{j) on the right-hand side contains all the 
terms acting on 'H±, including creation operators. The 
coefficients in (13) can again be found by taking commu- 
tators, for example: 



(14a) 



7j = [do,ad,out{j)] = ^i'iuv'^ $*d{j) 
Sj = [ai, ad,out(j)] = -^[fyj] vv^Tph- a*7j](14b) 

In general, a click recorded in the APD corresponds 
to the measurement of at least one photon. In the limit 
of low detection probability, the action of the detection 
is the subtraction of a single photon. Note that this as- 
sumption is practically always obeyed if the parameter j 
is in a continuum (like in the case of spectral filtering), 
since the probability to have two photons exactly in the 
same mode is then negligible. Henceforth we assume to 
be in this limit. In the Heisenberg picture, a photon 
detection then corresponds to the application of the op- 
erator a<j_out(i) to the initial state {i.e. to the vacuum 
state), followed by the normalization of the result: 

= PJ^'^aa^^MQ) (15) 
where Pj is the detection probability for the mode j: 

Pi = (0|aLut(i)ad,out(.7)|0) 

= {0\aUj)ad±{jm + \^j\' + \6j\' (16) 

= 0j;(j>i;'l' (/)d(j). 
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Below we will use that the vacimni expectation value of 
fidJ-U) = alij_{j)ad±{j) can be expressed as Pj - |7jP - 

Before continuing, it can be instructive to recall the 
concision allowed by the Heisenberg picture. In a single- 
mode problem, a photon-subtracted squeezed state is 
equivalent to a squeezed single-photon state: this case 
corresponds to the simple Bogoliubov transform Sout = 
cosh(r)ain -I- sinh.{r)al^, which directly gives a pure 1- 
photon state (after normalization) when applied to the 
vacuum. As the states do not evolve in the Heisenberg 
picture, they all can be considered as 'input' states; but 
when measured using output quadratures, this 1-photon 
state will appear to be squeezed. 

We arc going to use the same approach in the mul- 
timode case. One can first note that the conditioned 
state l^'(i)) in Eq. (15) is already a 1-photon state. This 
state, however, does not belong to 7^2 only, so that 
measurements output are not so obvious to compute. 
In fact, we are solely interested in expectation values 
{9{(ih,out, ajj out)) operators describing the output that 
is measured in the homodyne detector. Such expectation 
values can be written as 

{9{ah ,OUt5 out ) m (17a) 

= Tr{g{ah,out,&loutMj){i^j\}.{17h) 

In Eq. (17b), the trace can be separated into a trace over 
Ti.2 and a trace over TC±, and the latter does not act on 

the function 5(a/i,out) aj^ out)> whose expectation value is 
then: 

= T^2{9{ah,out,&l,JTr^\i^j){ijj\} (18a) 

= Tr2{5(afc,out,aI,out)Pj} (l^b) 
All quantities of interest can therefore be deduced from 
the input state reduced density matrix pj, acting in TI2, 
and the crucial advantage of mode reduction is to allow 
a simple expression for this matrix: Writing |0) = |00) ® 
|0)_L, where |00) and |0)_l are the ground states of H2 
and W_L, respectively, and using Eqs. (13,15), we directly 
obtain: 



Pj = Tr^|V,)(V',l 

= (1 -C,) |00)(00| +e,-4(j)|00)(00|a,(j), 
in terms of the modal purity 



(19) 



(20) 



and where 



with tan6'j = Sj/lj (21) 



is an operator that creates a single photon in a superpo- 
sition of mode and mode 1. The state (19) produced 
from a detection event in mode j is a mixed state, mix- 
ing vacuum and a single-photon state with weight ^j. 
Without conditioning or in the limit 0, we have 

Pj = mm- 



D. Wigner functions 

Squeezed vacuum. — Before determining the out- 
put Wigner function corresponding to the conditional 
state (19), it is instructive to first determine the Wigner 
function of the output state in the simplest experi- 
mental situation, where we ignore the photon detec- 
tor. The input state is then p = |00)(00|. We will 
make use of the standard Wigner functions of the vac- 
uum Wo{x,p) — exp(— r^)/7r and of single-photon states 
Wi{x,p) = (2r^ - 1) cxp(-r2)/7r, both with =x^+p'^. 
Clearly, W[p]{xo,po-, xi,pi) equals Wq{xq,po)Wo{xi,pi). 

In order to obtain the output Wigner function for the 
homodyne mode, we wish to express W[p] as a function of 
Xh.out, Ph.out (defined from the output homodyne mode 
fl/»,out given by Eq. (7)). This requires the introduction 
of another mode Si^out orthogonal to a/t,out) so that the 
transformation (ao,ai) (a/i.out, ai.out) is symplectic 
(i.e. commutation relations are preserved). Using the 
model of Fig. 2, one can choose Oi^out of the form: 



ai,out = P 



a ao 



+ vT+^ai. 



(22) 



This form is by no means unique, but this does not pose 

a problem since mode lout will eventually be integrated 
out. One can now invert the relations (7,22), thereby 
expressing a;o,i,Po,i as a function of x?(,out,P?t,out,a;i,out, 
and pi.out- After tracing over mode lout, which amounts 
to integrating over xi^out and pi,out5 we find the output 
signal entering the homodyne detector to be a squeezed 
vacuum state with a gaussian Wigner function 



Wo,sqz{x,p) 



1 



exp 



X 



(23) 



where x and p stand for x^^out and Ph,out and with vari- 
ances 



(24a) 
(24b) 



We can now invert Eq. (24) and rewrite the three mode- 
reduction parameters a, r], and (3 in Eq. (7) in terms of 
the variances, giving: 



V = 



a = 



2^V, + Vp + 2' 

2,/V.Vp - 1 
2^V^ + Vp + 2' 



(25a) 
(25b) 

(25c) 



In the following, we will keep writing a, 77, and (3 to 
shorten notation. It should be kept in mind, however, 
that Eq. (25) directly expresses these parameters in terms 
of the meastirable variances K,. p of the squeezed vacuum. 
In particular, (3 vanishes for minimal-uncertainty states. 
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Notice also that the parametrization for the mode- 
reduction parameters (25) is equivalent to the one in 
Eq. (12) in terms of squeezing parameters r and g. Thus 
r and g can be expressed in terms of the variances Vx,p, 
and vice versa. 

Photon- subtracted squeezed vacuum. — As for the 
squeezed vacuum, we now calculate the Wigner func- 
tion for the photon-subtracted squeezed vacuum, starting 
with the initial state (19). The mode (21) has a one- 
photon excitation in state (19). The orthogonal mode 
with creation operator aj_|_7r/2(i) "^^^ excited. Hence 
the Wigner function corresponding to the state (19) is 

W{xo,po;xi,pi) = (1 -Cj)W^o(a;o,Po)Wo(a;i,pi) (26) 

Note that in this expression, quadratures XQ. g.^^/2i 
Pej,ej+TT/2 can be easily expressed as functions of quadra- 
tures xo,i,po,i using (21). As before, the Wigner function 
for the output signal is found by using the symplectic 
transformation defined by Eqs. (7) and (22). By tracing 
again over the mode lout) we obtain (see appendix A) 

W,{x,p) = [c, + 2A~ + 2B,^ + D,^^ W^o.sqz. 

(27) 

The constants in this Wigner function are given by 

= Pf'hj{v + a) + Sjf3\\ (28a) 
Bj = P,-'hj{v-a)-6j(3f, (28b) 
Cj = l-Aj/V^-Bj/Vp, (28c) 
Dj = -8Pr^Im{^*Sj)rjp. (28d) 

Note that the Wigner function (27) of the photon- 
subtracted squeezed state differs from the Wigner func- 
tion of the squeezed vacuum Wo,sqz(a;,p) of Eq. (23) 
only because of the polynomial in x and p between the 
large brackets. The same quantities Vx.p as in Eq. (24) 
show up, with or without conditioning. Since in general 
W{x,p) > -TT-i [27], we find the condition Cj > -1. 

Averaged Wigner functions. — Practical detectors do 
not resolve with infinite precision when and where pho- 
tons are detected. We should therefore average over all 
possible microscopic states that agree with the detection 
record. We assumed in Sec. II C that the average num- 
ber of photons detected per pulse in the APD is much 
smaller than one. Averaging over unresolved detection 
events is then equivalent to averaging over single-photon 
subtraction events. 

The Wigner transformation of the density matrix is a 
linear transformation. Therefore, the averaged Wigner 
function W{x,p) is simply obtained by replacing Aj...Dj 
in (27) by A...D, with the notation 

X = P,-IJ2p,Xj. (29) 
j 

Here Ptot is the sum of the probabilities Pj of micro- 
scopic states that agree with the detection record. From 



Eq. (28) it follows that averaged quantities A...D involve 
sums like J2j l7iP) J2j l<^jP or J2j in the following 
we will write the respective averages as |7p, jjp or 7*5. 

The Wigner function (27) of the photon-subtracted 
state and its detection-averaged version have a very gen- 
eral significance. Before going to Sec. Ill, devoted to 
a practical implementation of these results, let us finish 
with some reflections on the detection modes. 

E. Detection modes 

All the previous results were derived through the use of 
a set of detection modes ad{j). The coefficients entering 
in the ave r aged Wigner function W involve quantities like 
-ftot, |7p, |<5p or 7*5, in which the detection operators 
only appear through their projection operator: 

n = Y.U3)$\{3)- (30) 

3 

Hence one would find the same predicted averages if one 
would employ a different set of detection modes that has 
the same associated projector. This projector describes 
how the setup filters the signal before it enters the photon 
detector. For example, a time-domain filtering system 
will be described by H = Ht, where Ht can be written 
in terms of a set of modes k{t) that are labeled by the 
time t: 

UT{t,t')=T{t)5{t-t'), (31) 

where T{t) = 1 when the APD is switched on, and 

T{t) = otherwise. To give another important exam- 
ple, a spectral slit can be described with a set of modes 
a{u)) with an associated projector 

Un{LJ,Lu')=T{Lj)S{uj~Lj'), (32) 

where T{uj) = if the frequency lu is filtered out, and 
T{u)) = 1 otherwise. 

Above we have assumed that the filtering of the sig- 
nal after its production in the DOPA was included in the 
Bogoliubov transform U . Below we will give an alterna- 
tive description, in which U is separated into the trans- 
formation due to the production of the squeezed light in 
the DOPA, and the subsequent filtering before detection. 
This alternative description will enable a more straight- 
forward comparison with the empirical model in Sec. HI. 

So, instead of the input-output transform Eq. (6) for 
the photon detection operator, we now write 

ad,out(j) = $li{j)uf{ua. + va*), (33) 

where the Bogoliubov matrix Uf accounts for filters (as 
the filters are passive, we have Vf = 0); this matrix u/ is 
of course unitary, even if the filters can present losses. In 
fact, losses will be modeled using beamsplitters, where 
the lost energy is reflected into auxiliary non-relevant 
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modes. These modes do not interact with the rest of 
the experiment (i.e. they are unaffected by BogoUubov 
transform u,v) and will not reach the APD. But all the 
other modes, referred to as relevant modes, should be 
considered as detection modes, and we then have: 



E^ci(i)<?d(i) = n, 



(34) 



where 11^ is the projector onto the subspace of relevant 
modes. Let be the projector on the non-relevant 
modes. As the latter are unaffected by the transform 
u,v, we have Ilr{ua, + va.*) = n^a. Inserting the relation 
Ilr + Ilr = 1 into (33) then leads to 

ad,out{j) = $dij)'^f^r{u a. + va*) + $l{j)uf'n.ra. (35) 

Here the last term on the right annihilates vacuum, and 
commutes with annihilation operators like ciq or oi: this 
term will add no contribution to the results of the pre- 
vious subsection. The only change therefore consists in 
the substitution (f UrU^^(f, so that the operator 11 in 
Eq. (30) should be replaced by 

n' = Uru\ J2 Mj)^liU)ufIlr = F^F, (36) 



where F = Uj-UfUr and where we have used standard 
properties of projection operators (11 = 11'^, 11^ = 11). 
The operator F represents the action of the filters re- 
stricted to the subspace of relevant modes. If F is a 
projector, such as IIt of Eq. (31) or Iln of Eq. (32), then 
the effect of 11' is the same as of 11 in Eq. (30). This can 
be easily understood: it is equivalent to say that the fil- 
tered modes are blocked, or that they are first redirected 
into auxiliary modes, and then blocked. 

The operator 11' of Eq. (36) is a more general quantity 
than n in Eq. (30), however, since 11' need not be a pro- 
jection operator. It can for instance account for partial 
absorption of the modes. In that case, the spectral trans- 
mission T{uj) in Eq. (32) can assmne any value between 
and 1, to account for filtering systems more complex 
than a simple spectral slit. 

Furthermore, the above expressions can simply be gen- 
eralized to situations where several filters are used. For 
example, if a spectral slit IIq is followed by a time-domain 
filter IIt, then the above expressions for F and 11' be- 
come F = IItIIo, leading to 11' = EnllTlIo. 



III. APPLICATION 

At this stage, we have a complete description of the 
final state starting from the Bogoliubov transform (1). 
The results of the previous section are generally valid, 
since we started with a multimode model that was left 
unspecified. Our purpose now is to establish a concrete 
link with the photon-subtraction experiment as described 
in Ref. [10], and to improve its analysis. 
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FIG. 3: Simplified experimental setup: a squeezed vacuum 
is generated by a DOPA, whore pliotons of frequency 2cj are 
transformed into pairs of photons of frequency lj (in the ex- 
periment, the central frequency wo corresponds to a central 
wavelength of about 850 nm). The output signal of the DOPA 
is sampled by a beam splitter with low reflectivity R. If a 
photon is detected by the APD, then ideally it has been sub- 
tracted from the squeezed vacuum. 



A. Photon subtraction experiment 

In the experiment by Ourjoumtsev et al. [10], pulses 
of squeezed light are produced. The setup is sketched 
in Figure 3. A squeezed vacuum, produced in a single- 
pass DOPA (a KNbOs crystal) by down-conversion of 
frequency- doubled femtosecond laser pulses, is sampled 
by a beam splitter with low reflectivity R = 1 — T . Two 
mode filters are placed in front of the APD: a spatial 
mode filter, that consists in a single-mode fiber, and 
a spectral slit of width Jl. If a photon is detected by 
the APD, then ideally it has been subtracted from the 
squeezed vacuum. This subtraction leads to a 1-photon 
squeezed state, which is very close to a 'Schrodinger kit- 
ten' state. Quantum state tomography with a balanced 
homodyne detector [27] allows the complete reconstruc- 
tion of this highly non-gaussian quantum state of light. 

The Wigner function (27) was derived assuming that 
the mode reduction was performed on the mode entering 
the homodyne detector. To relate our results to the em- 
pirical model discussed in the next subsection, we here 
choose to perform the mode reduction to the signal di- 
rectly after the DOPA. We model the DOPA using the 
scheme presented in Fig. 2, where the parameters r and g 
are linked to the Bogoliubov transform through Eqs. (9- 
12). For the calculation of the modes ad,out(.?) detected 
by the APD, the sampling beamsplitter and the mode 
filters can be separately added to this transform, as ex- 
plained in Sec. II E. As stated above, here we chose not to 
include the sampling beam splitter into the mode reduc- 
tion. The Wigner function (27) then describes the signal 
just after the DOPA. We therefore still need to account 
for this sampling beam splitter between the DOPA and 
the homodyne detection, as well as for other losses. For 
example, one usually accounts for imperfections of the 
homodyne detection by adding a fictitious beam splitter 
of transmission ?7hom just before the homodyne detection, 
where r^hom is tli(3 homodyne detection cfhcicncy. Both 
those beam splitters can be easily implemented by re- 
placing the variances according to: 



('") 



x/p 



1 = »7homr[14/p - 1]. 



(37) 
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and by multiplying A, B and D by ?7homT. Before going 
into detailed calculations, let us first recall the empiri- 
cal model that was proposed in Ref. [10] to account for 
experimental results. 



B. Empirical model 

It is \iscful to recall the empirical model proposed in 
Ref. [10] to explain the experiments, and to see by what 
assumptions our multimode model reduces to it. The 
DOPA is again modeled as in Fig. 2, producing the same 
squeezed vacuum. However, in the empirical model it is 
assumed that the detected photon is either in the homo- 
dyne mode with probability ^, or in an orthogonal mode 
with probability (1 — ^). In the latter case, the detection 
event is not correlated with the homodyne measurement, 
and one simply performs a homodyne measurement on 
squeezed vacuum. 

The output density matrix obtained with the empirical 
model is similar to our pj in Eq. (19). In fact, the two 
would be identical if the detected photon was only due 
to photons in the mode V'^ or from Ti.^. This is in gen- 
eral not the case, however, as there will be an admixture 
from Si, out in the photon detection operator. In fact, 
in order to completely account for the multimode nature 
of this experiment, the empirical model should be modi- 
fied in the way depicted on Fig. 4, with the insertion of a 
beamsplitter of amplitude reflection and transmission co- 
efficients p and T that allows interference between a;; out 
and ai^out- A photon detection event in such a setup can 
indeed be equivalent to the application of a\{j) to the 
initial vacuum [see Eqs. (19,21)], provided 



- cosh(r) = — , 

r tan Bj tan 



(38) 



with tan(^?o) = (3/a- These angles 6*0 and Oj are mixing 
angles that fix the probability amplitudes of detection 
of a photon of mode and of mode 1. Our angle 9j in 
general depends both on the squeezing properties of the 
light source, and on the filtering of the signal before the 
photon detector, whereas the empirical 6q only depends 
on the source. Within a narrow-filter approximation that 
will be detailed in next subsection, such a setup can also 
account for averaged quantities (29) with the use of an 
average angle 6 instead of 6j (see Eq. (54) in the follow- 
ing). 

This possibility of interference between the homodyne 
signal and the ai^out signal is the crucial difference be- 
tween the multimode and the empirical models: More 
interference makes the empirical model worse. The es- 
sential assumption of the empirical model is thus that 
the photon detection operator does not have a contribu- 
tion from ai.out- Then 7j and 5j could be replaced by a 
and (3, respectively, according to Eqs. (7) and (13), and 
the angle 6j in Eq. (19) by (in which case, according 
to Eq. (38), the beamsplitter BS{p,t) in the equivalent 
model in Fig. 4 can be removed). 
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FIG. 4: When one applies the mode reduction procedure to 
the multimode model for photon subtraction as sketched in 
Fig. 1, the resulting Wigner function is equivalent to the two- 
mode model depicted in the figure. The DOPA may in gen- 
eral be represented by an ideal single mode DOPA and a 
two-mode NDOPA as in Fig. 2, and the photon counting is a 
combination of dark counts as well as a coherent mixture of 
the two output ports of the NDOPA. Compared to the em- 
pirical model developed in Ref. [10], the only difference is the 
presence of the beam splitter BS{p^ t) which was not present 
in the empirical model. 



The empirical Wigner function can be easily deduced 
from Eq. (27) with the above replacements, and has 

the same form after the replacement of A - ■ ■ D by 



A ■ ■ ■ n 

^cmp ^cmp- 



Coefficients Aemp and Bemp are obtained 



by multiplying A, B by ^/^, and replacing 7 and 5 by a 
and /3, respectively. This gives 



-^emp — ^ 



-^emp C 



(39) 



1 A^jYip/^X 



The coefficient Ccmp is given by Cemp 
Beinp/Vp, and J^emp vauishes. 

The empirical model produces intuitive results. How- 
ever, it requires justification. If large spectral slits would 
be used, then the homodyne mode and many other or- 
thogonal modes would hardly be affected by the slit. If 
the detected photon could have come from many modes 
orthogonal to iph , then the modal purity ^ would be unac- 
ceptably low, and also a large admixture of ai^out would 
enter the detection signal. Indeed, some of us found ex- 
perimentally that the spectral slit should be as narrow 
as possible, while still allowing the detection of a sig- 
nal, in order to find the highest modal purities (see also 
Sec. IIIC). Consequently, narrow slits have been used in 
the photon-subtraction experiment [10]. Although it is 
obvious that filtering is necessary, the use of a narrow 
spectral slit before the photon detector does not make 
the empirical model automatically valid. A quantitative 
comparison of both models is therefore needed to test the 
validity of the empirical model, as given below. 
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C. Concrete multimode model 

Let us now develop a simple spatio-temporal multi- 
mode model for which the Bogoliubov transformation can 
be written explicitly. We assume that light propagation 
inside the DOPA is described by modes of the form 



A{f, t) exp(ia;ot — ik- r), 



(40) 



where the plane wave is exactly phase-matched, and 
where the amplitude A satisfies the slowly-varying en- 
velope approximation (SVEA). This approximation does 
not hold for all the light that exits the nonlinear crys- 
tal, but the homodyne mode is supposed to be phase- 
matched, and we will assume that the filters before the 
APD block the modes that are not phase-matched. We 
will furthermore neglect diffraction effects within the 
DOPA. In the basis {x,y,t), where x and y are spatial 
variables running on the DOPA's transverse plane, the 
u and V of the Bogoliubov transform (1) then become 
diagonal (see appendix B), and we have 

^■out{x,y,t) = u{x,y,t)a{x,y,t) + v{x,y,t)&' {x,y,t), 

(41) 

in terms of operators that we assume to have commuta- 
tion relations 

[a{x,y,t),a\x',y',t')]=d{x-x')S{y-y')d{t-t'). (42) 

Coefficients in the transformation (41) have the form 

u{x,y,t) = cosh[ql Ep{x,y,t)], (43a) 
v{x,y,t) = smh[qlEp{x,y,t)]. (43b) 

Here Ep{x,y,t) is the pump-beam amplitude, which we 

assume to be real-valued. The parameter q takes into 
accoimt the nonlinearity of the crystal and I is its length. 
One can allow for Group Velocity Mismatch (GVM) in 
the crystal by convoluting Ep by a rectangular unit gate 
of duration Tg, the time separation induced by the GVM 
after passing the crystal (see appendix B) . More precisely, 
this convolution should be made twice, as we also have 
GVM for the Second Harmonic Generation (SHG) of the 
pump beam. The u and v are real-vahied functions if the 
pump beam Ep is so, which is a valid assumption if there 
is no frequency chirp. The homodyne mode tph{x, y, t) = 
eh{x,y,t) will also be taken real- valued in the following. 

The homodyne signal is then given by a/i^out 
i e/iaout, where integration over x,y,t is implied. Mode 
reduction now starts with the identification 

r]ao = j dxdydteh{x,y,t)u{x,y,t)a,{x,y,t), (44) 

which is Eq. (8) specified for our spatiotemporal model. 
The mode-reduction parameters are now given by spatio- 
temporal integrals, for example 

,^ = /d.d„d<.J(.,„.«)„n.,,,«) (45a) 

•qa = j dxdydte\{x,y,t)u{x,y,t)v{x,y,t).{45b) 



Further parameters can be found analogously. 

Averaging over photon detection events. — We have 
seen in Sec. II E that all averaged quantities can be ob- 
tained through the determination of the operator H de- 
fined in (36) when the filters are separately added to the 
Bogoliubov transform. In the considered experiment two 
filters are used: a rectangular spectral slit, that can be 
described using (32); and a monomode fiber that selects 
a single spatial mode (j)s{x,y), which we suppose to be 
real-valued, and therefore corresponds to the projector 
(l)s{x,y)(ps{x' ,y'). Note that in this experiment detec- 
tion times are unknown at the scale of pulses duration, 
so that there is no time-domain filtering. (In the analy- 
sis one should average over all possible photon detection 
times.) We therefore have to use 

H(a;, y, w; x', y', J) = (j)s{x, y)^s{x', y')T{w)5{u) - w')> 

(46) 

where T{u)) = rjcR if a; enters into the spectral slit, oth- 
erwise. Here R is the sampling beamsplitter reflc^ctivity, 
and rjc accounts for all other losses in the conditioning 
arm (APD efficiency, optics losses ...). As can be seen 
from Eq. (40), the field amplitudes are defined around a 
central frequency ujq (or 2ujo for the DOPA pump beam, 
see appendix B), so that the frequency lu = for the 
amplitude Fourier transform corresponds in fact to this 
central frequency; in that way, a rectangular spectral slit 
well centered around this central frequency can be de- 
fined as w e [—17/2, i7/2]. The operator H defined in 
Eq. (46) should be applied to Fourier-transformed mode 
functions: 



(47) 



Using Eqs. (16), (46) and (47), the average total photon 
detection probability per pulse then becomes 



47r2 



dxdydt v'^{x,y, t)(pl {x,y), 



(48) 



where the average is taken over the detection modes 
4'd,j{x,y,t). The expression (48) increases linearly with 
the filter width fl, but is valid only if fl is small enough 
to warrant the SVEA. Furthermore, when conditioning 
upon a click in the detector in Sec. II C, we assumed that 
-Ptot «C 1, an assumption that can now be tested with the 
explicit formula (48). 

In general, the parameters 7 and S appear in W as 
product sums like 77* = ^ ITjPi or 7^* = J^j^j^j- 
our concrete model, the evaluation of detection-averaged 
coefficients in the Wigner function involves integrals of 
the type: 



47r2 



n/2 



"^"^ j dxdydx'dy' ' ■ 



a/2 



dui f{x,y,u))h* {x',y',uj) 



_ (49) 
For example, the average 77* is found by substituting 

both / and h in (49) by the time-domain Fourier trans- 
form of the function Chuv (ps{x,y,t). [Here and in the 
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following, wc abbreviate products like /(x, y, t)h{x, y) by 
fh{x, y, t).] Since / and h are Fourier transforms of real- 
valued functions, the corresponding integrals (49) are 
rcal-valucd as well. Hence, all mode-reduction param- 
eters and coefficients in the Wigner functions are also 
real-valued. In particular, the coefficient D in W van- 
ishes, see Eq. (28d) . There is no difficulty to numerically 
evaluate integrals (49) and we will do that below, but let 
us first focus on an additional approximation that can 
considerably simplify these results, without becoming in- 
accurate. 

Narrow-filter approximation. We previously dis- 
cussed the experimental observation that the spectral 
slit should be as narrow as possible. Another simpli- 
fication is possible in that case, that simply consists in 
neglecting in the integrals (49) the frequency dependence 
of the mode profiles within the narrow width fi, i.e. 
f{x,y,oj) f{x,y,0) for < Let us recall that 

this value w = corresponds to the central frequency luq 
of the pulses, where the real- valued amplitudes present a 
maximum. This has two consequences: first, the presence 
of this maximum justifies a zero-order Taylor approxima- 
tion, provided the spectral width CI is much smaller than 
the spectral width of the pulses (typically 2tt/t, where 
r is the pulse duration); second, if all functions involved 
in (49) present a maximum at w = 0, then the narrow- 
filter approximation generates an upper-bound for these 
integrals, and therefore for quantities like |7^|, |i5^| or 
the modal purity ^ (see 20,29). This approximation will 
be applied and tested in Sec. HID, dedicated to the nu- 
merical results. This approximation brings the following 
simplification in the integrals (49): 



dxdy fix, y, 0) j dx'dy' h*{x', y', 0) (50) 
dxdydtf{x,y,t) [ dx'dy'dt' h*{x' ,y' ,t'). 



rjcRn 
47r2 

rjcRn 
47r2 



Evidently, wc end up with separate integrals over / and 
h, and using the definitions (14) we obtain the averages 
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2-K T] 



J dxdydtchuv (f>s{x,y,t) (51a) 

In the narrow-filter approximation, averages of products 
are simply given by products of averages, 77* = |7p, and 
7(5* = 7^*, etc. Essentially in the f2 ^ limit the filter 
removes any temporal information about the time the 
photon was emitted from the DOPA. The photodetection 
then corresponds to a single mode with w = 0, regardless 
of the average over detection times. 

We therefore find for the photon-subtracted squeezed 
state an average Wigner function of the form (28), with 
coefficients 



B = P,-l\l{Ti-a)-m\ 



(52a) 
(52b) 



and C =1- A/Vx - B /Vp and D = 0. Using the same 
substitution in (20,21) we can also introduce the averaged 
modal purity 



and the average angle 9 defined by: 
tan^ = 5/7. 



(53) 



(54) 



Constant profiles. — Before dealing with a more re- 
alistic case, it is interesting to focus on the case of 
constant profiles. Let us assume a constant value for 
Bh, Ep, and (ps within a space-time support of volume 
S. The normalization of the homodyne mode implies 
= 1/S, and equation (43,45) leads to r] = u, a = v, 
and (3 = 0. As /3 vanishes, the mode a\ is no more de- 
fined, and Eqs. (14b, 51b) cannot be used anymore. In 
fact the mode-reduction procedure now leads to an effec- 
tive single-mode model rather than a two-mode model. 
The homodyne mode is now in the single-mode space 
Hi spanned by ao,aJ. One can simply put 5j = 
in (13), and hence 5 = 0. This leads to the average 
angle ^ = ^0 = 0. The modal purity becomes ^ = 1, as 
it should for a single-mode model. Most importantly, we 
find C = — 1 , which according to Eq. (28) corresponds to 
the most negative value for the Wigner function at the 
origin, W{0,0) = -I/tt. 

So, with constant profiles and a narrow filter slit, we 
recover from our multimode model the single-mode de- 
scription for photon subtraction experiments. Since this 
limit leads to the most negative Wigner function, it repre- 
sents the ideal limit for producing states for QIP applica- 
tions, at least according to our simple multimode model. 
This shows that the multimode nature essentially appears 
through the mode distortions due to the non-constant 
space and time profiles of the pulses. The central role of 
gain-induced distortions is particularly clear with regard 
to the multimode nature of the squeezed vacuum pro- 
duced by the DOPA: assuming a constant pump field, 
that is assuming no gain-induced distortions, is in fact 
enough to obtain (3 = 0. Let us now return to a more 
realistic model, taking into account these profiles. 

Gaussian profiles. — As a more realistic simplification, 
we assume gaussian profiles for the various fields. For 
instance, we write the homodyne field eh as 



eh{x, y, t) = Bhfi exp ( — — ) exp ( -2 



(55) 



where w is the beam waist, r the duration of the gaus- 
sian pulse, and e^.o a normalization constant. The pump 
beam is usually obtained by SHG, in a crystal pumped 
by a beam identical to the homodyne beam. In the low- 
est order of the SHG process, the profiles of Ep and 
have the same shapes, so we can assume another gaussian 
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profile: 



Ep{x,y,t) = EQep{x,y,t) 



= Eq exp 



exp 



where one expects the pump pulse duration to be tp = 
r/\/2. However, if GVM is taken into account, this gaus- 
sian profile (56) must be convoluted by rectangular gates. 
In practice, such convolutions lead to beam profiles that 
are still very close to gaussians. 

The final gaussian profile to be introduced here is the 
spatial mode 4's{x,y) of the filter in front of the APD. 

It is the LPqi mode of a monomodc fiber that is well 
approximated by a normalized gaussian of waist Wf. 



D. Numerical results 

Here our goal is twofold: first, to compare our multi- 
mode analysis with the empirical model that was used be- 
fore to analyze photon subtraction experiments. Second, 
by exploring our multi-parameter multimode model, we 
look for parameter regimes that are best suited for pro- 
ducing states with the most negative Wigner functions. 

Expansions in pump field. — For our numerical work it 
is convenient to write all fields as Taylor expansions in 
the pump field Ep. Prom Eq. (43) it follows directly that 



In the same way we have 
(56) r]a = ^CmPm-, 

m 

' -Plot = T}cR^ dmQmi 

m 

Vv^sr R 



u'^{x,y,t) 
uv{x,y,t) 
v'^{x,y,t) 



i + i cosh[2qlEp] = ^ 6„e^\ (57a) 



2 2 
1 
2 

2 ' 2 



^sinh[2gZ£p] = 5^c„e^, (57b) 

m 

-l + l cosh[2qlEp] = ^d„e?(57c) 



with Ep = Eoep{x,y,t) as in Eq. (56). This defines 
the constant coefficients bm,c„i, and dm- For qIEq < 1 
these expansions converge qiiitc quickly. We then only 
have to insert relations (57) into the various integrals for 
an efficient numerical evaluation. For instance, we can 
rewrite (45a) as: 



(58) 



with 



Prr 



da;dydi e^e™ = 



2x/2r-2 



WpTp 



(mu;2 + 2w%)^/mT^ + 2Tf,' 

(59) 



1] 



= —-\ 2^ d,nRr, 



/3 



(60a) 
(60b) 

(60c) 

(60d) 



with 

Qm 
Pm 



J dxdydt^le^ = 
^ j dxdydt(l)sepeh 

7r-3/4 



2V2m{mwf + 2wy) 

(61b) 



2\ ■ 



After fixing parameters, these expansions in the pump 
field can be readily used for numerical evaluations. 

Fixing basic parameters. — First we fix some parame- 
ters of our multimode model in order to present numerical 
results and to see how much our analysis differs from the 
one in Rcf. [10], where filtering before the photon detec- 
tion was not modeled explicitly. We take w = 1.2wp and 
a transmission T = 90% of the sampling beam splitter. 
Moreover, we fix W{ = w/1.5, which is compatible with 
the coupling efficiency into the filtering monomode fiber 
(approximately 80%, see Ref. [10]). 

Regarding efficiency of homodyne detection, the mode 
cih considered in Sec. 11 was defined as the mode that 
perfectly matches the local oscillator of the homodyne 
detection, in other words the matching efficiency equals 
unity by definition in our model. The transmission of the 
optics and the photodctcction efficiency together lead to 
an overall efficiency of homodyne detection 77hom- We put 
Vhom = 0.93, in agreement with Ref. [10]. 

As stated above, if GVM is taken into account, the 
almost gaussian profile of the pump pulse is convoluted 
twice by a rectangular gate with time window Tg. A 
KNbOs crystal of length I = 100 ^im has Tg = 120 fs. 
For an initial duration of the homodyne pulse t ~ 150 fs, 
the convolutions indeed lead to a nearly gaussian beam 
profile with Tp « r. We assume the identity rp = r in 
the following. 

Negative Wigner functions. — As stated in the Intro- 
duction, the global minimum of a Wigner function is the 
standard figure of merit for the nonclassicality and 'non- 
gaussianity' of the corresponding state. After subtrac- 
tion of a single photon, the Wigner function W{x,p) is 
always most negative in the origin (since D = 0). Fig- 
ure 5 shows how Ty(0, 0) depends on the squeezing factor 
s = exp(— 2r). The most negative values are obtained in 
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FIG. 5: (color online). Minimal value of the Wigner function, 
W{0,0), as a function of squeezing parameter s = exp(— 2r), 
which is varied by changing the quantity qlEo. The narrow- 
filter approximation was made for the spectral slit. Fixed 
parameters: pulse parameters w = 1.2wp, tp — t = 150 fs, 
transmission of sampling beam splitter T — 90%, efficiency of 
homodyne detection ryhom = 0.93. 
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Spectral filter transmission 

FIG. 6: (color online). Minimal value of the Wigner function, 
W{0, 0), as a function of the spectral slit transmission for the 
homodyne mode. Calculated for s = 0.56, using a complete 
evaluation of integrals (49), i.e. without the narrow- filter ap- 
proximation. The narrow- filter approximation is well satisfied 
for low transmissions. 



the low-squeezing limit s —> 1. This can be understood 
as there is less gain-induced distortions in that case. 

In Ref. [10], the best experimental results (highest 
modal purities) were obtained for s — 0.56. For this 
value of s, which can be selected by choosing the right 
value for the quantity qIEq, we obtain g = 0.50 and 
1^(0,0) = —0.034; the latter value is close to what was 
observed in [10], without correction for the detection ef- 
ficiency. At this stage, it can be interesting to compare 
this result, obtained using the narrow-filter approxima- 
tion, with a more accurate calculation based on a com- 
plete evaluation of integrals (49). Figure 6 presents the 
numerical results obtained for 1^(0, 0) at s = 0.56 as a 
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FIG. 7: (color online). Minimal value of the Wigner function, 
i.e. W{0,0), as a function of wp/w, where wp is the waist 
of the pump field and w is the waist of the homodyne field. 
qlEo is fixed such that squeezing parameter s = 0.56. Other 
parameter values as in Fig. 5. 



function of the spectral slit transmission for the homo- 
dyne mode. (This transmission can be increased by mak- 
ing the spectral slit width Q larger.) A minimal value is 
clearly reached for low transmissions, justifying a poste- 
riori the use of narrow spectral slits in the experiment 
of Ref. [10]. Since for low transmissions, W^(0,0) does 
not differ much from its minimal value, the narrow-filter 
approximation that we made in Sec. Ill C gives accurate 
results. 

Figure 7 predicts the behavior of W^(0, 0) when varying 
the size of the pump beam. Experimental values for the 
widths were related by w = 1.2ajp [10]. Fig. 7 clearly 
shows that one can await a high increase of the negativity 
from a larger pump beam. This result was intuitive, as 
there is less gain-induced distortions in that case, but is 
here quantified. This can motivate the use of amplified 
pulses [28], in order to have a spatially broader pump 
beam (i.e. with larger wp), but with the same intensity. 

Comparison with empirical model. — Measured nega- 
tive Wigner functions were interpreted in [10] using the 
empirical model as introduced in Sec. IIIB, where the 
photon is subtracted in the 'good' mode a^^out with prob- 
ability ^, and where the state is left in the squeezed vac- 
uum with probability (1 — ^). As explained before, the 
main difference between the empirical and our models 
is that the conditioned state in the former corresponds 
to the single-photon initial state Og^ with tan^o — P/ct, 
while it corresponds to Og in the latter, with tan^ = 5/7. 
These angles and 9 are mixing angles that fix the prob- 
ability amplitudes of detection of a photon of mode and 
of mode 1. Figure 8 shows 6*0 and ^ as a function of the 
squeezing parameter s. Clearly, and 9 do not differ 
too much, and by less than 10% for s = 0.56. 

Another important difference between our model and 
the empirical model is that modal purities in our model 
are fixed by the relation (20), whereas the parameter ^ 



12 




0.2 0.4 0.6 0.8 1 
s 



FIG. 8: (color online). Mixing angles 9 (solid line) and Oq 
(dashed line) of modes and 1, as a function of squeezing pa- 
rameter s, which is varied by changing qlEo- Other parameter 
values as in Fig. 5. 



in the empirical model is a free parameter. This freedom 
can be used to fit the data, i.e. to have Acmp = A or 
Bcmp = B. It is however not a priori possible to fit both 
parameters A and B (52a, 52b) from (39) using only the 
fitting parameter ^. In neither model should the vari- 
ances Vx and Vp be considered as free fitting parameters 
of the photon-subtraction experiment, at least their val- 
ues should agree with the values for V^^p obtained by 
homodyne measurements of the squeezed vacuum. 

In our model the mixing angles 9j and their average 9 
take into account the filtering of the signal that is used for 
conditioning. In the empirical model, the corresponding 
angle is independent of the filtering. Thus it is to 
be expected that this inaccuracy of the empirical model 
will lead to optimally fitted modal purities ^opt in the 
empirical model that are systematically lower than the 
average modal purity ^ in our model. This is indeed 
what we find for the curves in Figure 9: the best fit in 
the present example is obtained for ^opt = 0.87, a value 
that is indeed smaller but still close to ^ « 0.91. The 
high quality of this fit (with an error less than 1.2%) is 
directly linked to the fact that in the present case 9o « 9. 

There is a possibility to improve this result if g is con- 
sidered as a fitting parameter as well. We obtained an 
error of less than 0.4% between the Wigner functions for 
Copt — 0.89 and gopt = 0.53, i.e. for a value of g that dif- 
fers by 6% from the value given by the multimode model. 
In other words, if g has a great influence on Acmp and 
Bcmp in (39), it has a very low impact on Vx, Vp] the 
change from g = 0.50 to gopt = 0.53 modifies the val- 
ues oi Vx, Vp by only a few 10"'^, and for this reason it 
is very difficult to accurately measure g from squeezed 
vacuum [26]. These considerations explain why the em- 
pirical model can fit experimental data so successfully; 
even when 9q is not equal to 9, the parameter g gives a 
supplementary freedom for fitting. 
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FIG. 9: (color online). Sections along (a;,0) and (0,p) of the 
Wigner function W{x,p) (solid line), and corresponding best 
fits using the empirical model (dashed line) with ^opt = 0.87. 
Parameter values as in Fig. 5. Results of our model and the 
empirical fits almost overlap. 



IV. DISCUSSION AND CONCLUSIONS 

We have introduced a straightforward and physically 
intuitive procedure that we call 'mode reduction' to sim- 
plify the multimode description of squeezed light to the 
bare essentials. For photon-subtraction experiments, this 
means that the homodyne signal is reduced to an effective 
two-mode description and the detector signal requires one 
extra orthogonal effective mode. We derived the Wigner 
function of the homodyne signal conditional upon the 
detection of a single photon, and we also showed how to 
average over possible measurement outcomes. 

The general mode-reduction formalism was then ap- 
plied to a detailed model describing photon subtraction 
of gaussian spatiotemporal pulses of squeezed light. This 
model features many experimental parameters such as 
beam waists and duration of the pulses that can be in- 
dependently measured. Indeed, our model does not have 
free fitting parameters. This allows one to study in de- 
tail what are the crucial experimental parameters to pro- 
duce optimally negative Wigner functions with pulses of 
squeezed light. 

We compared our new model to the empirical model 
that was used before to analyze photon-subtraction ex- 
periments in [10]. In fact, the formulae for the output 
Wigner functions look similar. One crucial difference 
is that the empirical model does have a free parame- 
ter, namely the quantity called the modal purity. In our 
model modal purities also occur, be it with a slightly dif- 
ferent meaning, but they are fixed quantities. A good 
agreement between our model and experiments therefore 
gives more understanding than an accurate fit with the 
empirical model. 

We found that in the range of parameters of the mea- 
surements in [10], both our model and the empirical 
model are accurate. We reasoned that modal purities 
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in our model would be systematically higher, and in our 
numerical example we found this to be the case. The 
accuracy of the empirical model strongly depends on the 
availability of the free parameter. It was nevertheless a 
surprise in the theoretical analysis that the mixing angles 
9o and 0, describing the relative probability of measuring 
a photon in either one of two effective modes, differ at 
most 20% in a whole range of squeezing parameters. 

Our mode-reduction procedure is closely related to the 
analysis of photon-subtraction experiments of Refs. [19, 
20, 21]. One could express our mode-reduction param- 
eters in terms of elements of the covariance matrix of 
Refs. [19, 20, 21]. Our output Wigner function in Eq. (27) 
then reduces to the one in Ref. [21], but only in the 
special case that all our mode reduction parameters are 
real- valued so that Dj in Eq. (28d) vanishes. This we 
assumed for simplicity in Sec. IIIC. Our mode-reduction 
procedure is carried out in the Heisenberg picture. We 
think that our approach has some advantages. In our 
approach it becomes quite intuitive in what sense it goes 
beyond the empirical model of Ref. [10]. In our concrete 
multimode analysis, we include effects not considered in 
Ref. [21], such as the transverse beam profile, for which 
we found that wider pump beams lead to more negative 
Wigner functions. 

In conclusion, we presented a very concise model 
that can account for the multimode nature of projec- 
tive photon-counting measurements. It gives an intuitive 
picture of photon-subtraction experiments, close to the 
empirical model previously published. This multimode 
model therefore gives consistent results, in agreement 
with previously published experiments where pulses of 
light with negative Wigner functions were produced con- 
ditionally. Our model can be used to predict the changes 
in the output upon variation of experimentally relevant 
parameters, and to optimize the setup design. 
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Pi. out- It is however convenient to make a change of vari- 
ables so that the integral is over xi, pi instead of a;i^out) 
Pi.out- In this case the only transforms needed for this 
calculation is 



Xo 

Po 



rj — a 



l+/?2 

rj + a 



{Xh,out - Pxi) 
2 (P/i,out + PPi), 



as well as the transformation of the integral 

1 



/ 



da::i,outdpi 



I + /32 



da;idpi. 



(A3a) 

(A3b) 

(A4) 



One should then note that the Wigner function (26) is 
the product of a polynomial in x, p, and of a gaussian 
term exp (— i?^), with 

= xI+pI+xI+pI = Xg+pl+xl^^/2+Pe+w/2- (A5) 
With Eq. (A3), the exponent can be rewritten as 



Vx ( PXh,out\^ , ^lout 
Xl 77 1 + 



(?7 -I- a)^ 

[rj — a)2 



Vx 



Pi 



\ 2 2 
Pph ,out \ Ph,out 



Vr, 



Vr. 



The integral (A4) with Eq. (26) as its integrand can then 
be found by replacing in the integrand the squares x\ and 
pI by 



(A6a) 
(A6b) 



P\ 









2T4 




(?? - af 




'2Vp ' 



and by replacing the first-order terms according to 

liXh,ont 



Xi 



Pi 



liph,out 



(A7a) 
(A7b) 



APPENDIX B: SLOWLY VARYING ENVELOPE 
APPROXIMATION 



APPENDIX A: WIGNER FUNCTION 

In order to derive Eq. (27) from Eq. (26), one has first 
to invert Eqs. (7,22), leading to 

/- (^'A.oMt fiat oMt,^^^ 

ai = Vl + /32ai,out - /3aI,ouf (^2) 

One should then replace in Eq.(26) a;o,i,po,i by 
Xh,out,Ph,out,Xi,out, and pi,out, and integrate over Xi^out, 



The goal of this appendix is the derivation of the local 

Bogoliubov transformation (41). We assume that inside 

the DOPA the pump pulse with an angular frequency 

2u!o travels at a speed Ug,2wo > with negligible absorption. 

This field can therefore be written as 

z -* 
iEp{x,y,t 6t) exp{2iu!Qt — ik^uio • (^1) 

where 5t is an arbitrary time delay and where 'i' is a 
purely conventional phase factor. Let us write the probe 

beam as 

E{r, t) = A{r, t) exp(iwoi - «C ' 0> (B2) 
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where the phase-matching condition A:2cjo = '^k^,, is as- 
sumed to be satisfied. By using the SVEA in MaxweU's 
equations, neglecting diffraction terms and considering 
the first-order dispersion, we obtain 

r)A 1 r> A 7 

^ + %-=qEp{x,v,t-^-bt)A\ (B3) 

where A = A{f,t). The substitution of t — z/vg^^jo by t 
then leads to 

dA 

— (f, t) = qEp{x, y,t-Dz- 6t)A*{f, t), (B4) 

where D = v~l^^ — w^ig is the GVM. With the assump- 
tion that Ep is real-valued, the solution to Eq. (B4) be- 
comes 

Aout = cosh[ qlFp ] Ai„ + sinh[ qlFp ] A^^, (B5) 



where the (a;, y, t)-dcpendence was suppressed. The ef- 
fective pump field Fp is given by 



Fp = - dzEp^- dTEp{x,y,t-T), (B6) 

' Jo Tg Jst 



with Tg = Dl the time separation induced by the GVM 
after crossing the crystal. Eq. (B6) shows that the ef- 
fective pump field Fp is a convolution of the pump field 
Ep with a rectangular unit gate of duration Tg, which for 
5t = —Tg/2 is centered around the origin. As the quan- 
tized version of Eq. (B5), we then find Eq. (41) of the 
main text. The pump field Ep of the main text is to be 
understood as the effective pump field Fp derived here. 
Evidently, Fp Ep in the limit Tg ^ (no GVM). 
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